Dynamics, synchronization, and quantum phase transitions of two dissipative spins 
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We analyze the static and dynamical properties of two Ising-coupled quantum spins embedded 
in a common bosonic bath as an archetype of dissipative quantum mechanics. First, we elucidate 
the ground state phase diagram for an ohmic and a subohmic bath using a combination of bosonic 
numerical renormalization group (NRG), analytical techniques and intuitive arguments. Second, 
employing the time-dependent NRG we investigate the system's rich dynamical behavior arising from 
the complex interplay between spin-spin and spin-bath interactions. Interestingly, spin oscillations 
can synchronize due to the proximity of the common non-Markovian bath and the system displays 
highly entangled steady states for certain nonequilibrium initial preparations. We complement 
our non-perturbative numerical results by exact analytical solutions when available and provide 
quantitative limits on the applicability of the perturbative Bloch-Redfield approach at weak coupling. 
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I. INTRODUCTION 

A quantum system is never completely isolated from 
its environment which results in noticeable effects such as 
decoherence, dissipation and entanglement.^ One promi- 
nent example embodies a two- level (spin- 1/2) system in- 
teracting with a collection of harmonic oscillators, the 
so-called spin-boson model. ^""^ The latter displays a rich 
behavior ranging from damped Rabi oscillations to local- 
ization in one of the two states, and has been widely stud- 
ied as a paradigm of quantum dissipation and quantum- 
to-classical transitions.' As it constitutes the elementary 
unit of a quantum computer (qubit), much work was re- 
cently directed toward understanding and controlling the 
dissipative spin-boson dynamics in nonequilibrium situ- 
ations such as time-dependent external fields.*'"^" The 
model is of particular importance because it may be 
implemented in a variety of different experimental con- 
texts, for example, the tunneling of defects in solid-state 
systems, ^^ electron transfer in chemical reactions^^'^'' or 
qubit designs based on the Josephson effect. ^'^"^^ Other 
systems that are described by the spin-boson Hamil- 
tonian are trapped ions,^' quantum emitters coupled 
to surface plasmons,'^''* and the cold-atom quantum dot 



setup 
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Further variants of spin-boson models involve 



two-level atoms interacting with a single quantized mode 
of an electromagnetic cavity. ^"^"^^ 

The environmental influence on the phase coherence 
between the two spin states is of crucial importance in 
the field of quantum computing, as it sets a limit to the 
timescale where coherent quantum logical operations can 
be performed. In this context, it is essential to extend 
the system to multiple two- level systems (or qubits), as 
operations involving two-qubits, e.g., the CNOT gate, 
are required to obtain a complete set of quantum log- 
ical operations. In addition, the presence of a second 
spin allows to address the competition between spin-spin 



and spin-bath interactions and the resulting interplay be- 
tween quantum control and dissipation. 

In the present article, we investigate such a general- 
ization of the single spin-boson model and consider two 
quantum spins {cri, cr2} that are coupled to each other 
via an Ising-type coupling and interact with a common 
bath of harmonic oscillator modes, as described by the 
Hamiltonian (sec also Fig. 1) 
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We set the reduced Planck constant h = 1. Here, Ci'l'^ 
are the usual Pauli matrices describing the two spins and 
bk is the bosonic annihilation operator of the bath mode 
with frequency Wfc. The free spin part of the Hamiltonian 
contains the tunneling amplitudes Ai^2, bias fields ei.2 
and the bare Ising interaction constant K. The effects of 
the bosonic environment on the spins are fully captured 
by the bath spectral density^"''^ 

J{uj) = tt'^XIS{uj - ujk) = 2'n:aui''ujl~''6{ujc - uj)9{uj) , 
fe 

(2) 
which we assume to behave as a power law w'* (s > 0) up 
to the cutoff frequency Wc. Hereafter, we will be studying 
exponents in the range ^ < s < 1, where the case s = 1 
(s < 1) refers to an ohmic (subohmic) bosonic bath. The 
strength of the coupling to the bath is characterized by 
the dimensionless dissipation constant a > 0. 

For the spin-bath interaction, for simplicity, we use 
identical coupling constants A^ for both spins. This cor- 
responds to the case where the spins are spatially close to 
each other. Specifically, we assume their separation di2 
to be smaller than the shortest wavelength of the bath 
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FIG. 1: Two quantum spins-|, cri and CT2, coupled through 
an Ising interaction K. The spins are also entangled, via their 
cr^-components, to a common reservoir of bosonic oscillator 
modes with frequencies cuk- The bath is characterized by the 
spectral density J{oj) = 2-Rau>''ijj].'~"6(ijJc — u])9{lo), where s = 
1 (s < 1) refers to an ohmic (subohmic) bosonic environment. 



excitations: di2 "^ ^c — Vg/ojc, where Vs is the sound 
velocity in the bath.'^"'^^ 

There are several reasons for considering an Ising-like 
coupling ^(Tfa2 between the two spins. First, there 
are experimental situations where such an SU(2)-broken 
coupling is realized, for instance in capacitively coupled 
quantum dots where the operators a^ describe charge 
states on the dot.^^"''^ Other examples are the cold-atom 
quantum dot setting, trapped ions and superconducting 
qubits. Second, since the bath couples to the a^ com- 
ponent of the spins, it automatically induces an indirect 
(ferromagnetic) Ising interaction between the spins which 
is mediated by a coherent exchange of phonons. This re- 
sults in a renormalization of K to Kr — K ~ Aaud s. 
Therefore, even for zero K , the spins are Ising-coupled. 
We note that, in general, the bath induced interaction 
decays with the spatial distance between the spins di2 
on a lengthscale given by Ac ~ uj~^ .'^^''^'' 

The two-spin boson model allows to address the com- 
petition between spin-spin entanglement, characterized 
for instance by the concurrence, and spin-bath entan- 
glement, characterized for instance by the entanglement 
entropy. 
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The entanglement entropy also contains 
information about the coherence between different spin 
states.^ We will show below that for a particular initial 
preparation, the system exhibits a non-trivial steady- 
state, where the spins are strongly entangled with the 
bath while maintaining coherence between different spin 
configurations. 

Whereas for some experimental realizations the de- 
scription of independent bosonic reservoirs is appropri- 
ate, e.g., in the case of quantum dots coupled to indepen- 
dent leads, ^* there are others, where the spins couple to 
a common bath, e.g.., the cold-atom'^'^^^^ and trapped ion 
setup. ^^ Here, we assume a common bath because we are 
mostly interested in studying the competition between 
the coherent and dissipative parts of the interaction in- 
duced by the bath, leading to dynamical spin synchro- 
nization and highly entangled steady states. The other 
situation has been addressed for instance in Refs. 36-38. 



In the following, we aim to investigate not only the 
static properties of the ground state but also the nonequi- 
librium dynamics of the system, both for an ohmic and 
a subohmic boson bath. In the subohmic case, we 
mainly consider the experimentally relevant situation of 
s — 1/2.'^^''^'^ We apply the powerful non-perturbative nu- 
merical renormalization group (NRG).''^^'''' To solve for 
the dynamics of the system, we employ the recently de- 
veloped time-dependent NRG (TD-NRG),'*^''^^ that we 
compare to exact solutions, available at special points in 
the parameter space, and to the Bloch-Redfield master 
equation approach'*'^ at weak dissipation. 

The paper is outlined as follows. In Sec. II, we calcu- 
late the zero temperature phase diagram as a function 
of dissipation strength a and Ising coupling K, both for 
s — 1/2 and s = 1. As a reminiscence of the single spin- 
boson model, it contains a delocalized phase ((erf 2) = 0) 
for small dissipation and a localized phase ((erf 2) 7^ for 
£1^2 = 0"*") for large dissipation. We give a physically in- 
tuitive explanation for the asymmetry between the ferro- 
magnetic {K < 0) and antiferromagnetic (K > 0) regions 
of the phase diagram. 

In Sec. Ill, we investigate the critical properties at the 
phase transitions such as the behavior of the entangle- 
ment entropy across the transition, or the scaling of spin 
expectation values, that occurs for a subohmic bath. 

In Sec. IV, we explore the nonequilibrium dynamics 
of the two spins after a quantum quench of parameters. 
We typically polarize the spins initially by applying large 
bias fields along the z or a;-direction that we switch off at 
time t — 0. We begin our analysis in Sec. IV A with the 
exactly solvable case of zero transverse fields Ai^2 = 0, 
where we show that our TD-NRG results perfectly agree 
with the exact analytical solution. In Sec. IV B we inves- 
tigate the regime of weak spin-bath coupling, and com- 
pare TD-NRG to the commonly employed perturbative 
Bloch-Redfield approach. We give quantitative limits on 
the applicability of the Redfield method. In Sec. IV C, 
we find that, interestingly, the bath is able to synchro- 
nize spin oscillations via a coherent exchange of phonons, 
even at weak spin-bath coupling. This phenomenon is 
not captured in the Bloch-Redfield master equation ap- 
proach, where the backaction of the bath on the spins is 
neglected. This method thus fails to correctly describe 
the spin dynamics even in the perturbative regime. In 
Sec. IV D, we investigate the spin dynamics for vanish- 
ing (renormalized) Ising interaction Kr = and high- 
light similarities and differences to the single spin-boson 
model. We elaborate on the case of weak-dissipation in 
Sec. IV D 1, where we compute the quality factor of the 
damped oscillations. In Sec. IV D 2, we discuss the dy- 
namics at the generalized Toulouse point a = 1/2. In 
Sec. IV E, we examine the crossover to the regime of 
strong spin-bath coupling for general Ising coupling, and 
point out differences between the case of an ohmic and a 
subohmic bath. In Sec. IV F, we describe that a highly 
entangled steady state can emerge from the dynamics if 
the system is prepared far from equilibrium. We finally 
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FIG. 2: Phase diagram of the ohmic two-spin boson model 
as a function of dissipation strength a and Ising couphng K. 
Different curves correspond to different values of tunneling 
amplitudes Ai = A2 = A. For infinitesimal bias fields ei,2 = 



-fO~ 



the ground state of the system in the localized region 



is given by Itt) ® |^), where |fi) is a shifted bath vacuum [see 
Eq. (6)]. The dashed line indicates where the renormalized 
Ising interaction vanishes: Kr = 0. 



conclude in Sec. V, and leave the details of some of our 
calculations to the Appendix. 



FIG. 3; Phase diagram of the subohmic two-spin boson model 
with s = 1/2 versus a and K, and for different values of A. 
The dashed line indicates where Kr = 0. 



First, in Sec. II A, we present the numerically obtained 
phase diagrams. Then, in Sec. II B we perform a strong- 
coupling analysis that will provide us with a qualitative 
understanding of the underlying physics. 



A. NRG phase diagrams 



II. GROUND STATE PHASES 



In this Section, we employ the bosonic NRG^ 



to 



calculate the ground state phase diagram corresponding 
to the Hamiltonian in Eq. (1) as a function of dissipation 
strength a and Ising coupling K. We present the results 
for the ohmic, s = 1, as well as the subohmic case of 
s = 1/2. We point out similarities and differences to the 
situation of the single spin-boson model and to that of a 
two-spin model with two separate baths. 

Throughout this study we use the following parameters 
for our NRG calculations (we use the common notation) : 
a discretization parameter of A = 1.4, a total of N^fi = 
599 bosonic modes in the first iteration and Ni,^M = 6 in 
the following ones, while keeping TVlov — 200 low-energy 
levels in each NRG iteration. 

We obtain a qualitative understanding of the phase 
diagram by using the fact that the fast bath modes fol- 
low the spin dynamics adiabatically in the sense known 
from the famous Born-Oppenheimer approximation.^''^ 
The spins are dressed by the bath phonons, and as a re- 
sult the energy separation of the two lowest-energy spin 
states becomes renormalized. This situation is reminis- 
cent of the single spin-boson model. There, the tunnel- 
ing splitting A also becomes renormalized by the bath, 
and in the ohmic case, one finds a renormalized value of 
Ar = A(^)"/'^~"^ for a < 1 and a complete quench of 
the tunneling for a > 1, where the system is thus local- 
ized."''' 



Using the NRG, we have determined the phase diagram 
of the two-spin boson model in Eq. (1). We present re- 
sults for an ohmic bath'*'^'"'^ in Fig. 2 and for a subohmic 
bath with s = 1/2 in Fig. 3. Different curves correspond 
to different values of A/ujc- Here, we assume equal tun- 
neling amplitudes of the two spins Ai = A2 = A. Intro- 
ducing slightly asymmetric tunneling elements Ai 7^ A2, 
however, does not affect the location of the phase bound- 
ary much. Hereafter, we use units of the bath cutoff 
frequency, i.e., we set Wc = 1, and we shall be mainly in- 
terested in the case where both Ai^2 ^ ^c and £1,2 ^ Wc- 

As shown in Figs. 2 and 3, the two-spin boson model 
exhibits two ground state phases: a delocalized phase, 
where the spin expectation values {af 2) vanish in the 
ground state for £1^2 — >■ 0, and a localized phase, where the 
spins develop a finite magnetization (af) — (cr|) = ±m 
{m > 0) for infinitesimal bias fields ei_2 — 0^. Like in 
the single spin-boson model, the system is delocalized 
for weak dissipation and enters a localized phase upon 
increasing a. The phase boundary, however, now explic- 
itly depends on the Ising interaction constant K. 

Let us first focus on the ohmic model in Fig. 2. For 
ferromagnetic K < 0, the phase boundary only weakly 
depends on K and is located at ttc ~ 0.15-f C'(^), which 
is a much smaller value than in the single spin case, where 
the transition occurs at a^'^^ie = 1 _)_ 0(A) 2,3 p^j. anti- 
ferromagnetic AT > 0, we find that the delocalized region 
extends up to larger values of a and we observe that the 
phase boundary occurs at the line K = AauJc/s for larger 



values of K. At this value of K the renormalized Ising in- 
teraction Kr, which takes into account the bath induced 
ferromagnetic spin-spin interaction (—AaiOc/s), vanishes. 
We defer the derivation of this formula until Sec. II B. 

Let us now turn to the subohmic case in Fig. 3. It 
shows the same qualitative features as the ohmic one, 
however, the system enters the localized phase for even 
smaller values of a. On the ferromagnetic side K < 0, 
our results suggest that ac ~ + C'(^), in agreement 
with the single spin case.'^" For antiferromagnetic i^ > 0, 
the system again remains delocalized up to larger values 
of a and the phase transition occurs close to the line 
Kj. = 0. Note that K^ depends on the bath exponent s. 

We distinguish the two phases by applying small bias 
fields ei_2 = lO^^ojc and measure (o'f2)- The latter 
vanishes in the delocalized region, but remains nonzero 
(erf) = (erf) = —m [m > 0) in the localized part of 
the phase diagram. We have also applied an antiferro- 
magnetic bias field configuration ei = —62 = 10~^uJc to 
test whether the system can also localize in an antiferro- 
magnetic spin configuration {|ti)) |4-)t)}- Interestingly, 
however, we observe in Fig. 4 that the spins always lo- 
calize in one of the ferromagnetic spin states (Itt)) III)}- 
The system does not localize in any of the antiferromag- 
netic spin configurations. We provide a physical explana- 
tion for this phenomenon in Sec. II B. Results for (0-^2) 
as a function of a for both bias field configurations and 
different values of K are shown in Fig. 4. We observe 
that (17^2) remains zero up to a larger value of a (for 
fixed K) simply because the antiferromagnetic bias fields 
£1 = — e2 do not lift the degeneracy of the two ground 
states {|tt)i 14-4-)}- The location of the phase boundary 
does of course not depend on the infinitesimal fields. 
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FIG. 4: (0-^2) as a function of a for various values of K and 
A — 0.025 (jjc. Different bias field configurations are shown in 
the upper part (ferromagnetic, ei = £2 = 10~*a;c) and lower 
part (antiferromagnetic, ei — —62 = W~^u!c) of the figure. 
This plot shows that spins are always aligned in the localized 
phase. The expectation values (0-12) remain zero up to larger 
values of a simply because the antiferromagnetic bias field 
configuration does not lift the degeneracy of the ground states 
(Itt)) \ii)} iu the localized phase. 



of if to 



Kr=K- 



Aaujr 



(3) 



This is most easily derived by applying the polaron uni- 
tary transformation [/ — exp[— i(crf-|-cr|)^j. ^{b\ — bk)] 
to the Hamiltonian in Eq. (1), which yields for H = 



B. Qualitative understanding of the phase diagram 

From the previous considerations, immediately the 
questions arise why the phase diagram is not symmet- 
ric under the combined transformation of {K — >■ —K, 
(T2 —7- — C2}, and why the system cannot localize in one 
of the antiferromagnetic spin states {|ti), lit)}- 

In order to answer these questions, we perform a 
strong-coupling analysis which relies on the fact that the 
fast modes of the bath (wj, ^ A) adiabatically renor- 
malize the energy separation of different spin states." In 
physical terms, assuming that the bath oscillators fol- 
low the time evolution of the spins immediately (Born- 
Oppenheimer approximation), the spins are dressed by 
phonons with frequencies larger than A. Thus, transi- 
tions between different spin states are suppressed if they 
involve a readjustment of the bath excitations. We will 
consider the ferromagnetic and antiferromagnetic cases 
separately. 

Let us first note, however, that the bath induces a fer- 
romagnetic interaction between the spins, which renor- 
malizes the value of the Ising constant from its bare value 



-A 



H = Y.[f{.le^''+h.c.) + 



Kr 



J = l 



Y.'^kblbk, 



(4) 



k>0 



where the hermitian bath displacement operator reads 






(5) 



This form of the Hamiltonian makes explicit the bath 
induced ferromagnetic Ising interaction. In particular, if 
the bare Ising coupling is antiferromagnetic K > 0, the 
effective interaction changes sign at a dissipation strength 
of a = 1^. For larger values oi K > ojc, the phase 
transition occurs close to this critical value of a, as shown 
in Figs. 2 and 3. From H, we can also learn immediately 
that a spin flip is associated with a complex excitation 



of the bosonic bath into a coherent state IfJ) 



jn 



|0), 



where |0) is the ground state of the free bath part of the 
Hamiltonian Hb = '^k^kbkbk- 



With this in mind, let us begin our strong-couphng 
analysis of the phase diagram with the ferromagnetic sit- 
uation K < 0, and assume that \K\ 3> Ai 2 and zero bias 
ei_2 = 0. For Ai,2 = 0, the two lowest energy spin states 
are given by the two ferromagnetic states {jft); lii)}- 
If we now turn on the tunneling Ai = A2 = A, wc find 
that the energy splitting between the two lowest states is 
of the order 



2A2 

SE^—{n\-n), 



(6) 



where the coherent state \fl) = e |0) is also referred to 
as the displaced oscillator bath state. It occurs when all 
oscillators equilibrate in contact with spins that are held 
fixed in position |^^). In terms of the spectral density, 
the bath renormalized energy splitting becomes 



SE ^ 



2A2 
W\ 



exp 



duj 



JH 



pSE 



(7) 



where p ^ 1. To be consistent with the adiabatic renor- 
malization scheme, the energy splitting 6E shows up as 
an infrared cutoff for the oscillator frequencies that are 
summed over. Since the bath renormalizes the energy 
splitting to smaller values 6E < "r^, one can solve 

Eq. (7) iteratively.^'* In the case that SE is renormal- 
ized to zero, the ground state is doubly degenerate and 
the system localized. This situation, where the displaced 
bath states jfi) and | — J7) are orthogonal to each other, is 
known as orthogonality catastrophe."' If SE is renormal- 
ized to a nonzero value, the ground state is unique and 
the system delocalized. 

For a subohmic spectral density, the iteration pro- 
cess yields SE = for any positive value of a, and 
the system is localized as soon as a > 0. In the 
ohmic case, on the other hand, we find that as long as 
a < 1/2, the energy splitting renormalizes to the finite 

,2a/(l-2a) _^_ ,^ _ oa2, 



value SE = SEo{'^y'^"^ '"' where SE^ ^ 2A^/\K\ 
For a > 1/2, however, one finds SE = and the system 
is localized. The phase transition occurs at the critical 
value ac — 1/2. The same value was recently found using 
a variational treatment."'^ Let us remark that in the case 
of the single spin-boson model, one has to calculate the 
overlap integral (^| — ^) — exp[— ^ J^^ duj-^], which 
^'■^ This also implies that the delo- 



leads to al'^^^ie = 1/2 



calized phase in the two-spin case is characterized by a 
distinct Kondo scale compared to the single spin-boson 
model." 

Our NRG calculation, which goes beyond this simple 
approximation and the variational approach of Ref. 49, 
indeed shows that the critical value of a in the ferromag- 
netic regime only weakly depends on K. In the ohmic 
case, wc observe, however, that ac rather converges to 
ac{s = 1) « 0.15 for large \K\ and -^ — >■ instead of the 
approximated value ac = 1/2. In the subohmic case, on 
the other hand, NRG agrees with the predicted value of 
ac = as we find ads = 1/2) w for — — )• 0. 



We now turn to the antiferromagnetic situation K > 
0. Since we want to investigate the antiferromagnetic 
regime, we thus have to assume that Kr > {or K -^ 00 
for any value of a). Then, the two lowest energy states 
for zero tunneling (A = 0) are degenerate in energy and 
given by {|ti): lit)}- If wc turn on tunneling, the two 
states hybridize and the energy difference between the 
two lowest energy states reads 



2A2 A^ 



(8) 



where | ) is the unshifted bath vacuum. Hence, any 
nonzero value of A leads to a unique ground state, be- 
cause the quenching of the tunneling amplitude due to 
the bath does not occur for a total spin zero state [com- 
pare with Eq. (6)]. [This can also be interpreted as the 
disappearance of Kondo-type entanglement for a spin 
zero state."'] As a result, the system is always delocal- 
ized for an antiferromagnetic Ising coupling Kr > 0, and 
the phase transition to the localized state is shifted to 
much larger values of a necessary to compensate the an- 
tiferromagnetic spin-spin coupling constant K. 



III. PHASE TRANSITIONS AND SCALING 

In this Section, we investigate the behavior of the sys- 
tem close to the localization phase transition in more de- 
tail. It is known, that the transition is in the Kosterlitz- 
Thouless universality class for the ohmic system, '''"'^ but 
it is of continuous type in the subohmic case.'^'''''^ Since 
recent studies show that NRG is not well-suited to de- 
scribe the system correctly close to the transition for 
s < 1/2, ■''"'■''■'' we restrict ourselves to s > 1/2. 

In Sec. Ill A, we first study the behavior of the entan- 
glement entropy in the ohmic and subohmic system. We 
then examine in Sec. Ill B the scaling of the spin expec- 
tation values (erf 2) close to the phase transition in the 
subohmic system. We derive mean-field scaling relations 
for the critical exponents from an effective spin action 
functional, and compare the resulting exponents to the 
critical exponents that we extract from NRG. 



A. Static entanglement entropy 

The entanglement entropy £ quantifies the degree of 
entanglement between the spins and the bath. It is de- 
fined as^^ 



£ = -Ti-[ps log2 ps] 



(9) 



where ps ~ Tr^p is the reduced density matrix of the two 
spins. Here, Tr^ denotes taking the trace over the bath 
degrees of freedom, and p is the full density matrix of the 
spin-boson system. One finds that < £ < log2 4 = 2, 
where £ = in the absence of entanglement between spin 
and bath. In Fig. 5 we show results for the entanglement 




0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 






K/ui^ = 




f\ :' KK ?: , -0.5 — - 




\ ■ / \i ' / \ -0.3 




/ -1 i 0' / \ -■■ '■■ ■■ -0.1 -■«-■ 




/ / /;/ V i \ ^ ''■■ ■' • 0-0 "■ 




/ / f ? / i \ / '. ' 0.1 -•- 




\i :\ \ '■ ': : \ 0.3 -<y- - 




' '\ / ' ' X ''■■' \ 0.5 ■■.■■ 


/ 


si y • * \ .-'A i 




,. / '\ V \ .-' Q 




/ \ '• '■■ \ * 












// V\ i. K''\ "■. 


*: 


' A ''-• y\ \ 








B-' \..-\ ".. "X V^ a.. 


- 4 

&' 





0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 



FIG. 5: Entanglement entropy f as a function of dissipation a 
in the ohmic two-spin boson model, shown for different values 
of the Ising coupling K and Ai,2 = O.liOc-- The rapid drop to 
zero around Oc ~ 0.25 signifies the transition to the localized 
phase. The plateau for smaller dissipation indicates the loss of 
phase coherence at a « ac/2 similar to the single spin-boson 
case. The inset shows larger values of K where the (incoher- 
ent) plateau shrinks to a peak-like structure, indicating that 
coherence is lost only right at the phase transition. 



FIG. 6; Entanglement entropy f as a function of dissipation 
a for the subohmic two-spin boson model with s — 1/2. Dif- 
ferent curves are for different values of the Ising coupling K, 
and Ai,2 = O.lujc- The entropy £ reaches a maximum at the 
phase transition (see also Fig. -3), and falls off continuously to 
both sides of the transition. 



a nonzero expectation value {a^) 7^ (see Ref. 50 for the 
single and Sec. IV E for the two-spin boson model). 



entropy in the ohmic system as a function of dissipation a 
for different values of Ising coupling K. Like in the case of 
the single spin-boson model, the entanglement entropy is 
nonzero only in the delocalized phase and rapidly falls to 
zero at the phase transition. It reaches a plateau for a ~ 
ac/2, indicating that coherence is lost already before the 
system becomes localized. The plateau characterizes a 
region of maximal decoherence, where the spin dynamics 
is incoherent. This coherent-to- incoherent crossover is 
known from the single spin system, '^'■'''''''''' where it occurs 
exactly at the Toulouse point a = 1/2. In Sec. IV D 2, we 
discuss the equivalent of the Toulouse point in the two- 
spin model where it is located at a = 1/2 and K = 2uJc- 

Surprisingly, as we show in the inset of Fig. 5, the 
plateau shrinks considerably if we go to larger positive 
values of K > ujc- The plateau more and more resembles 
a peak-like structure. This indicates that the localization 
phase transition occurs much closer to the regime, where 
spin oscillations are coherent. Coherence is lost only right 
at the transition (similar to the subohmic case discussed 
below). This is different from the single spin case, where 
the incoherent regime extends between 1/2 < a < 1 and 
is thus much larger. 

Finally, we show in Fig. 6 that for a subohmic bath, the 
entanglement entropy rather reaches a maximum (peak) 
right at the localization quantum phase transition. This 
behavior is known from the single spin-boson system. ^ It 
signifies that the coherence of the spin oscillations (con- 
tinuously) decreases toward the phase transition. There 
is no region where the spin transitions are completely in- 
coherent. In fact, coherent spin oscillations of cr^(i) even 
persist into the localized phase, where they occur around 



B. Scaling of magnetization for subohmic bath 

In this Section, we investigate the scaling of the spin 
expectation values {af 2) (magnetization) at the phase 
transition in the subohmic system. For the single spin- 
boson system, it is known that the phase transition is 
continuous for s < 1, and scaling exponents have been 
extracted using NRG^'^^''"'^ and Quantum Monte-Carlo 
calculations.^^ Recently, it was realized that NRG is not 
well-suited to describe scaling correctly for s < 1/2 in the 
single spin-boson model. '''^ Therefore, we only consider 
exponents in the range 1/2 < s < 1. 

We proceed in the following manner. First, in 
Sec. IIIB 1, we derive an effective spin action functional 
by integrating over the bosonic degrees of freedom. From 
this action we determine, in Sec. IIIB 2, the scaling di- 
mension of the spin operators in a mean-field approxi- 
mation from which follow scaling laws. We compare the 
resulting mean-field values for the critical exponents to 
those that we have extracted from the NRG calculations, 
and find good agreement between most of them. On the 
one hand this justifies our mean-field approximation, but 
on the other hand it also shows that the NRG analysis 
goes beyond this approximation. 



1. Effective spin action functional 

An effective action functional S'off for the spins can 
be obtained by integrating over the bosonic degrees of 
freedom using a functional integral description.*^'' This 



can be done exactly, because the Hamiltonian in Eq. (1) 
is quadratic in bosonic operators. 

We start with the action of the full system S — Ss + 
Sb + Ssb, where Ss - /q" dr^^'^J^aJ(r) + ^a|(T)] + 
^crf{T)a2{T) depends on spin variables only, and Sb = 

/o ^'''Sfc ^fc(''')[g7 + ^k]bk{T) denotes the action of the 
free bath. The spin-bath interaction is described by 
SsB = Uo dTj:,j:'^^,X,a^^iT)[bliT) + bkir)]. Here, 
/3 = 1/T (fcs = 1) is an inverse temperature, r is an 
imaginary time variable and 6(t) are the usual complex 
boson coherent state variables. Note that in the end we 
will take the zero temperature limit which is well-defined 
in this formalism. ^"^ 

Integrating over the (complex) bosonic variables'^'^ 
/P[6*(r),6fc(r)]exp[-S'B - Ssb] = exp[-S"], leads to 
an effective spin action 5cff = Ss + 5". In the zero tem- 
perature limit, it takes the form 



2. Scaling analysis: comparison between mean-field and 
NRG exponents 

Below, we derive mean-field critical exponents from 
the effective spin action in Eq. (10), which we compare 
with exponents that we have extracted from our NRG 
calculations. 

To proceed, we resort to a mcan-field-like decoupling 
of the Ising term: ^^^af (r)a|(T) « ^[af{ai) + {cTf)al]. 
This term then acts as a single-spin detuning, depending 
on the expectation value of the other spin magnetization. 

Scaling of both spins will thus be identical and 
we can follow the analysis for the single spin-boson 
model. ^'■''^''''■^'®^ There, one employs the quantum-to- 
classical mapping of the spin-boson model to the one- 
dimensional classical Ising modeP'''^'^'^ 



-ff, 



classical 



/ , '^ij Si Sj 



m 



short-range 



(12) 



S'off = / '^'^V^]^^']^'^) + 'T^^'ji'^) + ~r"^i('^)'^2('r)f with long-range interaction Jij = J/\i~i\^+''. Here, 
"'0 j=i 5f — ±1 are classical Ising spins. There is an additional 
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drdr' 
167r 



2 

I do; J(c.)e--l--^'l {^ [a|(r) - a|(r')] }■ 

(10) 



i=i 



The effect of the bosons on the spins is twofold: first, 
the Ising interaction constant gets renormalized to Kr = 
K — AaoJc/ s by the term that is local in imaginary time. 
Second, the bath introduces dissipation as described by 
the last term in Eq. (10), which is purely non-local in 
imaginary time. Integrating over frequency w, we observe 
that this last term describes a long-range interaction in 
imaginary time 



Sf — ±.1 are classical Ising spins. There is an additional 
generic short-range interaction .ffghort-rango arising from 
the transverse field, but it is believed to be irrelevant 
for the critical behavior.'"'^''''''®'' The scaling dimensions of 
al 2 a-re thus solely determined by the dissipative term, 
and we find from the condition that the total action is 
dimensionless [S] = 1 that 



Ka] = T- 



(13) 



Here, we have used units of energy (or temperature): 
[t] = T^^ . From this follows the scaling dimension of 
the detuning and Ising constant as 



[e]=T- 



[K] = T" 



(14) 



dijjj{uj) 



o\t-t\ 



27raa;2r(l + s) 
(l+c^ck-T'|)i+- 



(11) 



where T[x) is the Gamma function and we have used 
an exponential cutoff for the spectral density J{uj) — 
2Traiul^'^u!^ exp[— cj/cJc] for convenience. 

Note that although the dissipative part still contains 
a term that couples the two different spins (at different 
times), this corresponds to a retarded Ising interaction 
and can thus be neglected compared to the equal-time 
contribution, if one is interested in ground-state proper- 
ties. More specifically, the retarded term is of the form 
J drdr' °'| _ Tn+s , which under a Fourier transformation 
becomes ^^ Iw^l^crf (a;„)cr|(— Wn). Thus, if we pass to 
real frequencies a;„ — > a; + i6 and take the low-frequency 
limit a; — >■ these terms can be neglected compared to 
the static Ising interaction part '^crj^(T)CT|(r). This rea- 
soning can also be justified by noting that one arrives 
at the same formula for the renormalized Ising constant 
by applying the polaron unitary transformation to the 
Hamiltonian in Eq. (1) as we have presented in Sec. II B 
(see Eq. (4)). 



In order to derive scaling relations, we need to make an 
ansatz for the impurity part of the free energy. Since the 
fixed point is "interacting" for s > 1/2,'''^'^^ we use 

Fi,„p = T/(|A - Ae|T-l/^ 6^-^ \K - K,\T--) . (15) 

This ansatz can be applied for s < 1 since the transition 
is continuous. Further, for a Gaussian fixed point, which 
occurs at s < 1/2, the reduced free energy would also 
depend on dangerously irrelevant variables. 

In this ansatz we have used that in a quantum phase 
transition, which occurs at T = 0, the distance to criti- 
cality is measured by the parameter deviation from the 
critical value of the most relevant perturbation, in this 
case I A — Ac|. Analogous to a classical system, where 
the correlation length diverges as a function of this dis- 
tance, here the correlation length in imaginary time 
obeys .J ~ |A — Ad""^ with the correlation length ex- 
ponent V. The dynamic critical exponent is formally set 
equal to z = 1 in this -I- 1-dimensional system. This 
defines a characteristic energy scale 



rji-M 



(16) 



Exponent 


s = k 


«-| 


s-^. 


6 


4 
3 


10 

7 


40 
19 


c 

Cmf 


0.5 

1/2 


0.2 
1/6 ~ 0.17 


0.1 
1/18 ~ 0.06 


/3 

/3MF(i^ = 1/s) 


0.5 

1/2 
1/4 


0.2 

1/6 

1/4^/2 ~ 0.18 


0.09 

1/18 

1/4^5 ~ 0.11 



TABLE I: Comparison of critical exponents as predicted by 
our mean-field analysis {5mf,Cmp, /3mf} and as extracted 
from NRG {(5,C,/3}. 



above which critical behavior is observed. ^^ 

Using the ansatz for the free energy given in Eq. (15), 
we can immediately infer from [eT~''] — [\K — Kc\T~'^] — 
1 that b = ii^ and k ~ s. If we define the critical 
exponents describing the scaling of the magnetization as 



H2> 
{<2) 






(17) 
(18) 
(19) 



we can derive mean-field scaling relations. For instance 
from Eqs. (13), (14) it immediately follows that 



Smf — 



1 



1 



c 



1 



MF 



2s 



We have to invoke Eq. (16) to arrive at 
Pmf = V 



(20) 



(21) 



If we use the result that v = 1/s for small s, derived in 
Ref. 57, we find that (^mf = Pm f ■ Close to s = 1 it is 
more appropriate to use 1/v — \J2{1 — s) as obtained in 
Ref. G7. The resulting values for the critical exponents 
are shown in Table I. 

Let us now compare these mean-field predictions of the 
critical exponents to our NRG results. Numerically, we 
investigate the cases s = {5, f, ij)}- After carefully de- 
termining the position of the phase transition, we keep 
all but one parameter fixed at their critical values, and 
study the scaling of the magnetization as a function of 
this remaining parameter. Typically, we find power law 
scaling over more than two orders of magnitude, and we 
find the exponents from simply fitting the slope in a log- 
log plot. We have checked that the extracted value of 
the exponent is independent of the position in the phase 
diagram where we cross the phase boundary. As an ex- 
ample, in Fig. 7, we show the scaling of {af) as a func- 
tion oi \K — Kc\- Different curves are for different values 
of the transverse field A, and we extract the value of 
({s = i) = 0.5, which is is perfect agreement with the 
mean-field prediction of Cmf{s = ^) = 1/2- In Ta- 
ble I we show a full comparison of the critical exponents 



fdK) <x \K, - A-|C< 

Ci = 0.48310± 0.00370 
0.48786± 0.00373 
0.49035± 0.00417 
0.49292± 0.00489 
0.49556± 0.00586 
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FIG. 7: Scaling of magnetization at the phase transition in 
the subohmic system with s = |. We fit (af) as a function of 
the Ising interaction against a power-law (ai) ^ \K — Kc]'' , 
and find ( = 0.5. Different lines represent fits using fi oc 
\K — A'd""' . Results of the fit and error bars for (^i, as well as 
the different values of Ai,2 (in units of cjc) used, are shown 
in the plot. 



derived in the mean-field approximation and extracted 
from NRG. Agreement is good for the exponents ( and 
/3 for all values of s (using the different expansions of z/ 
as a function of s). For the exponent 5, however, the 
agreement is not so good in the cases s = {|, j^}. Note 
that the trend is captured correctly and that S diverges 
as s — >■ 1 which makes it increasingly hard to extract its 
value numerically. 



IV. NONEQUILIBRIUM SPIN DYNAMICS 

Let us now turn to the dissipative nonequilibrium dy- 
namics of the two-spin boson model of Eq. (1). We will 
concentrate on the ohmic (s = 1) as well as the subohmic 
case of s = 2- To access the system's rich dynamical 
behavior arising from the interplay of spin-spin and spin- 
bath interactions, we employ the time-dependent numeri- 
cal renormalization group technique (TD-NRG), recently 
introduced by Anders and Schiller.'*^ Using this exten- 
sion of the standard bosonic numerical renormalization 
group, ''■^"^'' we are able to calculate the real-time evo- 
lution of an impurity observable as a reaction to a sin- 
gle sudden change of parameters. Since this method is 
non-perturbative as well as non-Markovian, it is capable 
to accurately describe the spin dynamics over the whole 
range of parameter values, including strong coupling. We 
note that a first dynamical study of the ohmic system in 
a very limited region of parameter space using analytical 
methods was given in Ref. 68. 

As common to all applications of the NRG to bosonic 
quantum impurity models, we have to restrict the maxi- 
mal number of bosonic degrees of freedom that are added 
in each step of the iterative diagonalization procedure 



performed within the NRG method. We have checked 
that this cutoff does not aher our resuhs. We use the 
same NRG parameters as for the equihbrium calcula- 
tions: a discretization parameter of A = 1.4, a total 
of Nbfl = 599 bosonic modes in the first iteration and 
Nb,N = 6 in the following ones, while keeping A^lcv = 200 
low-energy levels in each NRG iteration. For the TD- 
NRG calculations we have additionally averaged the real- 
time data using N^ = 8 independent NRG runs (z-trick 
averaging). For more details about the method, we refer 
the reader to Refs. 46,69. 

In the following, we discuss a number of different 
nonequilibrium situations. 

In Sec. IV A, we show that TD-NRG results perfectly 
agree with the exact solution that is available for zero 
transverse field Ai_2 = 0, where the Hamiltonian only 
contains the z-component of the spin operators. 

In Sec. IV B, we focus on the case of weak spin-bath 
coupling and compare TD-NRG to the commonly used 
perturbative Bloch-Redfield method. We provide quan- 
titative limits at which dissipation strength this method 
begins to fail. 

We discuss, in Sec. IV C, the fascinating phenomenon 
of dynamical synchronization of the spin oscillations in- 
duced by the bath. Most importantly, this feature oc- 
curs even at weak spin-bath coupling and synchroniza- 
tion can thus be observed over many oscillation periods. 
It relies on the coherent exchange of bath excitations be- 
tween the two spins, which gives rise to the bath induced 
part of the Ising interaction. The phenomenon cannot 
be observed within the Bloch-Redfield master equation 
approach, where the backaction of the bath on the spins 
is neglected. 

In Sec. IV D, we investigate the spin dynamics for van- 
ishing (renormalized) Ising coupling Kr = 0. Qualita- 
tively, the system behaves like a single spin-boson model 
for < a < 1/2, where it exhibits damped coherent os- 
cillations. The quality factor of the oscillations, however, 
is smaller in the two-spin case as the damping is stronger. 
Yet most importantly, for larger values of a we find that 
the two spins remain delocalized for Kj. = up to a dis- 
sipation strength as large as a = 1.5 in the ohmic case. 
The single spin-boson model, in contrast, becomes local- 
ized at a = 1, where the spin remains frozen in its initial 
state. In Sec. IV Dl, we first elaborate on the region 
< a < 1/2, and use an approximation that is known 
to be equivalent to the Non-Interacting Blip Approxima- 
tion (NIBA)^'^" in the single spin case. It allows us to 
understand the dynamics qualitatively. In Sec. IV D 2, 
we then focus on the generalized Toulouse point a = 1/2 
and Kr = 0, where erf 2(i) decays purely exponentially. 
We show that one obtains slightly different decay rates 
for the single and two-spin boson model. We qualita- 
tively explain this difference by employing a bosonization 
mapping to a fcrmionic resonant level model. In the sin- 
gle spin case, the fcrmionic model can be solved exactly. 
For two spins, however, the fcrmionic model contains an 
additional interaction term that stems from the Jordan- 




FIG. 8: Comparison of {o-f 2(i)) between exact solution in 
Eq. (23) (dashed) and TD-kRG result (solid) for different 
values of a and bath exponents s = 1 (upper part), s — ^ 
(lower part). Parameters used are ei,2 ~ Ai,2 = 0, K — 0, 
LUc = 1, and a as specified in the plot. 



Wigner transformation of the spins and impedes an exact 
solution. 

We discuss the spin dynamics at large spin bath cou- 
pling in Sec. IV E. Comparing the ohmic and subohmic 
cases, we find that while coherence is lost prior to local- 
ization in the ohmic system, the spins exhibit oscillations 
even inside the localized regime for a subohmic bath, a 
feature only recently discovered'''^ in the single spin-boson 
system. 

Finally, as presented in Sec. IV F, an interesting situa- 
tion arises if we prepare the spins in an antiferromagnetic 
initial state at a location in the phase diagram which 
corresponds to a localized (ferromagnetic) ground state. 
Following the spin's dynamics over time, we observe a 
non-trivial steady-state, where the spins are highly en- 
tangled with the bath while developing and maintaining 
coherence between the two antiferromagnetic spin states. 
We give a simple physical explanation for this behavior. 



A. Decoherence without transverse field 

In this Section we discuss a specific case where one can 
exactly solve for the (non-trivial) dissipative spin dynam- 
ics of the two-spin boson model. We compare the exact 
solution to the TD-NRG results and find perfect agree- 
ment, which provides another validation of this powerful 
method in the strong coupling regime. 

For vanishing transverse fields, A 1,2 = 0, the Hamilto- 
nian in Eq. (1) takes the form 



■i Z 

i/[Ai,2 = 0] = ^{^[6,+^Afc(fe: 



fc>0 



h)]} 



K 



'1"2 



(22) 



fc>0 
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which only contains the z-component of the spin opera- 
tors. Thus, the spin dynamics is non-trivial only if the 
initial state of the spins contains a transverse component 
(in the x or y-direction). For instance, spins that are ini- 
tially polarized along the x-direction, {(7^{t = 0)) = — 1 
for j = 1,2, will undergo damped oscillations in (aj) as 
exactly described by'^ 



K(0) 



K, 



-Y 





rQifi)! 




r-Q2(i)i 


cos 


TT 


exp 


TT 



(23) 

with functions Qi{t) = J^ dco J (to) sin urt and Q2{t) = 
J^ dujJ{u})[l — cosuit]. For an ohmic spectral density 
with exponential cutoff, they read^ 

'1 , 



Qi(t) = 27ratan ^ w^i 
Q2{t) = TTa\n[l + uj^t^] . 



(24) 
(25) 



Note that Eq. (23) can also be derived using the polaron 
transformation of Sec. II B. In Fig. 8 we compare our 
TD-NRG results with this exact analytical prediction for 
bath exponents s = ^ and s = 1, and we find perfect 
agreement between them. 
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t (in units of uj^^) 



FIG. 9; Comparison of the results for {crf,2(t)) from TD-NRG 
(solid) and Bloch-Redfield approach (dashed) for K — ei,2 = 
0, Ai,2 = O.luJc in the perturbative regime aln ^ <^ 1. Up- 
per (lower) part shows the case s = 1 {s — |). Deviations be- 
tween the two solutions are visible aheady for q In ^ = 0.01. 
Note the beatings in oscillations due to the bath induced Ising 
coupling. 



B. Breakdown of Bloch-Redfield description 

In this Section, we compare the results from TD- 
NRG with those from the commonly employed Bloch- 
Redfield'''''^''*^ formalism at weak spin-bath coupling. We 
give quantitative limits on the applicability of this per- 
turbative and Markovian technique. 

Let us briefly outline the Bloch-Redfield approach to 
dissipative spin dynamics. The time-evolution of the spin 
reduced density matrix ps — Tr^ (p) , where Tr^ denotes 
tracing out the bath degrees of freedom and p is the full 
density matrix of the spin-boson system, is given by the 
Bloch-Redfield equations 

PS,ab{t) = -i^abPS,ab{i) - ^ Rabkl P S ,kli^) ' (^6) 

k,l 

Here, a,b,k,l G {1, . . . , 4} label the four eigenstates (with 
eigenenergy Ea) of the free spin part of the Hamiltonian 

Hs = E'=i [^^j + T^|] + f '^1^1 ' and ojab = Ea-E, 
are transition frequencies. For zero bias ei.2 = 0, the 
eigenenergies are given by E12 = =F^-/2 and E^^i — 
=Frj+/2 with n± = v^(Ai± A2)2 + X2/4. The relevant 
transition frequencies for which {a\al 2I&) 7^ 0, read 0^41 — 
UJ23 — f2 and 1^42 = wis ~ S with il ~ ^{^+ + ^-) and 
6 = 5(^1+ — Jl_). For nonzero bias £1^2 7^ 0, one can 
easily diagonalize Hs numerically for a specific choice of 
parameters. 

The Redfield tensor Rabki describes the effect of the 
bath onto the spin dynamics in the Born-Markov approx- 
imation. ^^ The real part of Rabki describes the damping 
induced by the bath, and the imaginary part the renor- 
malization of the transition frequencies, up to second or- 
der in the spin-bath coupling constants {Afe}. 



The Redfield tensor is explicitly given by golden rule 
transition rates and reads 



Rabki = Sbl 22 ^arrk + ^o.k X, ^ Irrb ^ ^ Ibak ^ ^ 



(-) 
Ibak 
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The golden rule rates at temperature T = 1/ /3 are calcu- 
lated to 

FJJj, = / ■" [coth{l3huj^j/2) T 1] 



Ibak A 

i^lbak 



V / duj 



0.2 - 4. 



[coi\\{l3fujj/2)ujij ^ uj] , 



where uiij = uib for the plus rates F;, A and uj, 



.(-) 



^i] 



(28) 



Wafe 



for the minus rates F^jj^j.. Here, we have defined the 
transition matrix element 



^Ibak 



^IM^l.ak 



^2,lb^2,ak 



(29) 

and a spectral density that is antisymmetri- 
cally continued to negative frequencies J{^) = 
s\g'i\{ijj)'iTa\ujYijj]r'^9{ujc — |w|). At zero temperature, the 
real part of the rates becomes 



ReF 



(±) 

Ibak 



A 



Ibak 



Jiuj.j) , 



(30) 



where again ojy = Uak for the plus rate and Wy = ujib 
for the minus rate. Note that Eq. (30) vanishes unless 
LOij > 0. The principal part integral in the imaginary 
part of the rates can be performed analytically, and for 
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FIG. 10: Comparison of results for ((Tf_2(i)) obtained with 
TD-NRG (solid) and Bloch-Redfield approach (dashed) for 
K = 0.2 Wc, £1,2 = 0. Other parameters are as in Fig. 9. 
Deviations between the two solutions are visible already for 
a = 0.005. 



instance in the ohmic case and for a Drudc bath cutoff 

J{uj) = 2Traoj/{l + ^), we obtain 



imr 
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Ibak 
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Ibak naUJijU!^ 
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(31) 



In all our calculations, we use the corresponding expres- 
sion for a hard bath cutoff [see Eq. (2)] which turns out 
to look more complicated than Eq. (31), but leads to the 
same results as long as ujc is the largest energy scale in 
the system. 

In the Redfield approach, the total density matrix is 
assumed to always factorize into a spin and a bath part. 
Further, by taking the long time limit in Eq. (28), any 
reversible energy exchange between spins and bath, and 
thus any back action of the bath on the spins is ne- 
glected. Therefore, the Redfield approach does not cap- 
ture the bath induced non-dissipative spin-spin interac- 
tion —AauJc/s = Kr — K correctly. As we will discuss 
in Sec. IV C this has the important consequence that the 
phenomenon of a bath induced dynamical synchroniza- 
tion of the spin oscillations cannot be observed within 
the Redfield approach. 

In this Section, we want to focus on a symmetric 
setup of the two-spin system (Ai = A2 = A) and 
zero bias £12 = 0. To determine the breakdown of the 
Bloch-Redfield description, we compare in Fig. 9 our TD- 
NRG results for (erf 2(^)) with Bloch-Redfield solutions of 
Eq. (26) for an ohmic and a subohmic bath. Both results 
agree for very weak spin-bath coupling aln^ < 0.01. 
However, already at aln^ = 0.01 we find significant 
differences. They are more pronounced in the subohmic 
case, and grow with the coupling strength. Even in the 
absence of a direct Ising coupling term K = 0, we ob- 
serve beatings in the oscillations due to the bath induced 



FIG. 11: Synchronization of two spins with different spin-flip 
terms A2 = |Ai — 10~^a;c by weak coupling to a common 
bath. There is no direct coupling between the spins K — 0. 
The upper part of the figure shows the uncoupled case a — 0. 
The middle part shows the ohmic case and the lower part 
the subohmic one with s = |. We use the same strength 
a = 8 ■ 10~* for the different bath dispersions, which lies in 



the perturbative regime: a In 



A2 



lO-^' < 1. 



Ising interactions Kj.. In Fig. 10 we show results for a 
system with a direct Ising coupling, where the beatings 
are stronger. Here, wc find significant differences between 
the TD-NRG and Bloch-Redfield results to occur already 



for a — 0.005 or a In 



Ai 



0.01. 



In summary, since the Redfield approach does not cor- 
rectly account for the bath induced Ising interaction, its 
breakdown occurs not just when a In ^^^ « 1, but already 



for 



Here Wy is a (nonzero) transition fre- 



quency of the system which is of the order {fi, d}. Since 
LUc ^ Wjj the breakdown of the master equation descrip- 
tion occurs for much smaller values of a compared to the 
single spin case, where it takes place when a In -^^^ sa 1. 



C. Synchronization of spin dynamics 

In this Section, we address how the coupling of spins to 
a common bath can be employed to obtain a dynamical 
synchronization of spin oscillations. Notably, this fea- 
ture occurs already at weak spin-bath coupling, where 
the bath induced decoherence is small. It provides an 
alternative technique to synchronize the dynamics of a 
two-spin system, when a strong direct coupling of the 
spins is unavailable. 

Let us start with two free and uncoupled spins {K — 
a ~ 0) that are driven by different tunneling amplitudes, 
say Ai = 2A2. As shown in the upper part of Fig. 11, 
the spins will then undergo undamped Rabi oscillations 
with frequencies Ai and A2, respectively. We now con- 
sider a weak coupling to the bath in the perturbative 



12 



regime where a In -^^ <^ 1 . The frequency corrections 
in Ai 2 are smah in this case. However, the bath induced 
Ising interaction Kr — —Aauic/ s can still be compara- 
ble to Ai^2, because it scales with the (large) bath cutoff 
frequency Wc- In this case, where Kj. and Ai^2 are of 
the same order of magnitude, the bath is capable of syn- 
chronizing the spin oscillations as depicted in the two 
lower parts of Fig. 11 for an ohniic (middle) and a sub- 
ohmic bath with s — 1/2 (bottom). The synchronization 
is more complete for the subohmic system, since there 
is an increased number of slow oscillator modes present 
and the induced Ising interaction, which scales as s~^, is 
twice as large (for the same value of a) . 

The two oscillation frequencies {^2,(5} that occur in 
Fig. 11 can be calculated from the free spin dynamics 
of Eq. (1), if we set the Ising interaction K equal to its 
rcnormalized value K = Kj. — —Aaud s. 

For zero bias ei_2 = 0, the free spin part in Eq. (1) 

^crfal with eigenvalues 
= =F^+/2, where fi± = 
■\/(Ai ± A2)^ -I- K'^/4:. We can find the spin dynamics 
from (erf 2(*)) — Tr5[/95(i)CTf 2], where, in the absence 
of the bath, the spin density matrix psit) evolves in 
time according to the von-Neumann equation of motion 
PS = —i[Hs,ps]- With initial condition ps{0) = |-^i)(ii|, 
we find for j = 1, 2: 



reads Hs = Y.'j^i ^^J 
Ei^2 = T^-/2 and £^3,4 



{^])^J2ps^-b{oy 



'{b\a^\a) 



2 A'^^' cos nt + 2 A'g'hos St , 



iU). 




0.002 0.004 0.006 

K (in units of ujc) 

1,2) 



0.008 



0.01 



FIG. 12: Oscillation amplitudes yl^ 'g as a function of Ising 
coupling K (see Eqs. (37) and (38)). Other parameters read 
Ai = 2A2 = 2 ■ 10~'^ujc- The two spin expectation values 
{<^i,2(i)) evolve according to Eq. (32). 



i[Hs,(jf 2(^)1 in Laplace space. ''^'^ One finds that 
A(if + Ai + A2) 



(^I(A)) = 



(A2 + f72)(A2+^2) 

A(f + A2 + A2) 



(35) 



(36) 



(A2 + r!2)(A2 + ,52)' 

which yields Eq. (32) in real space. We identify the 
amplitudes A^^ 'g ' as the respective residues of Eqs. (35) 
and (36) at fl and 6. Explicitly, they read 



(32) 



(1) 


±[-i^2^4(^2_ 


-Al)]+w 


(37) 
(38) 


(2) _ 
A ~ 


Aw 


-Al)]+w 

1 



where a, 6 label the eigenstates of Hs and tUab = Ea~ Eh 
are the transition frequencies. They obey 0^41 = 0123 = ^ 
and UJ42 — '^is — S. The two oscillation frequencies that 
appear in Fig. 11 are thus given hy fl = ^(^+ -l-ri_) and 
6 = 5(ri_|_ — r2_). For the other transitions, we find that 
the matrix elements (6|cr||a) are equal to zero. The two 
oscillation amplitudes are given by (j = 1, 2) 



where w = ^\K'^ + 4(A2 + A2)]2 - 64A2A2 and the up- 

fi 2) 
per sign relates to Aq_ ' . Synchronization sets in when 



1(1 



1(1 



A^^ ~ Ag which occurs for an Ising interaction strength 



of 



K 



A? -All 



(39) 



(i) 
n 

U) _ 



Ps,4i(0)(lk||4) -fp5,23(0)(3|a||2) (33) 

Ps,42(0)(2|af|4)+ps,i3(0)(3|a||l). (34) 



They are shown in Fig. 12 as a function of Ising cou- 
pling K, and are responsible for the synchronization phe- 
nomenon. At K — 0, the first spin oscillates with fre- 
quency fl{K = 0) = Ai and Aq = 1. The second spin 

(2) 
oscillates with frequency 6{K = 0) = A2 and A^ — 1. 

As we increase K, the amplitude Ag increases, while 

Aq decreases, and A^g remain almost the same. For 
large K :^ Ai_2, both spins oscillate with frequency 

(5~2AiA2/X. ' 

In fact, we can derive analytical expressions for A^ 'g if 
we solve the Heisenberg equation of motion for af 2 {t) = 



The damping rates of the oscillation amplitudes are 
proportional to J{V,) and J{S), respectively [see Eq. (30)]. 
Since J(w) ~ aw** they are small for the parameters in 
Fig. 11, and the synchronized oscillations can be seen 
over many periods. 

We like to emphasize that the synchronization effect 
cannot be seen in the Bloch-Redfield master equation 
treatment, which does not correctly account for the bath 
induced Ising interaction Kr (see Sec. IV B). The damp- 
ing rates in Fig. 11, however, agree with the ones calcu- 
lated with the perturbative Redfield approach. 

In summary, the bath induced Ising interaction scales 
with the bath cutoff frequency like K^ ^ auJc whereas the 
bath induced damping is proportional to afl'' and aJ*, 
where {f2, 6} are spin transition frequencies. A common 
bath can thus synchronize spin oscillations at weak cou- 
pling if the bath cutoff frequency is large: ujc 3> 51, S. 
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D. Vanishing Ising interaction Kr = 0: similarities 
and differences with the single spin-boson model 

In this Section, we investigate the spin dynamics along 
the Une Kj. = in the phase diagram, i.e. for a vanish- 
ing renormahzed Ising interaction (see Figs. 2 and 3). At 
first sight, one might expect that the dynamics would be 
identical to that of two uncoupled spin-boson systems. 
However, as can be seen from the perturbative treat- 
ment in Sec. IV B already, the two spins do not decouple 
from each other even to linear order in a, and the golden 
rule rates in Eq. (28) contain the terms al hi(T2 ak ^^'^ 
<J2ib^iak- Quahtatively though, as displayed in Figs. 13 
and 14, the spin dynamics of the single and two-spin bo- 
son models agree for Kj. — and small a. 

In the ohmic case of Fig. 13, we observe a crossover 
from damped coherent oscillations at < a < 1/2 to 
incoherent behavior at a > 1/2 in both models. In the 
two-spin case, however, we find stronger damping due to 
the terms proportional to crftT| mentioned above. This 
results in a smaller quality factor of the oscillations. We 
compare the quality factor of oscillations for the single 
and two-spin boson system computed with TD-NRG and 
Bloch-Redfield in Fig. 15. A detailed discussion of the 
dynamics at the special Toulouse point a — 1/2 is given 
separately in Sec. IV D 2. 

If we further increase a, we surprisingly observe that 
the two-spin boson model does not enter the localized 
phase (for Kr = 0). Unlike the single spin case, the two 
spins remain delocalized up to values of a > 1. Our time- 
dependent numerical results in Fig. 13 show that {erf 2} 
relax to zero even for values as large as a = 1.5. We note 
that this is in agreement with the NRG phase diagram in 
Fig. 2, which shows that the position of the localization 
phase transition converges toward the line Kr — from 
the side where Kr < 0. 

In Fig. 14, we show the same comparison between sin- 
gle and two-spin boson model for a subohmic bath with 
s = 1/2. As before, for increasing dissipation the coher- 
ence of the spin oscillations is lost more rapidly in time, 
and a comparison of the quality factors of the single and 
two-spin boson system is presented in Fig. 15. Again, 
the system does not localize along Kr — 0, which is in 
agreement with the NRG phase diagram of Fig. 3. 

In the following, we first derive in Sec. IV D 1 a decou- 
pling approximation that is equivalent to the well-known 
Non-Interacting Blip Approximation NIBA for the single 
spin-boson dynamics. In this approximation the spins de- 
couple completely for Kr — 0. It allows us to understand 
the spin dynamics along the line Kr = qualitatively. 
In contrast to the case of the single spin-boson model, 
however, the approximation does not give quantitatively 
correct results for two spins, not even to linear order in a. 
The reason is that dissipative second-order processes that 
involve both spins [see Eqs. (29) and (30)] are neglected. 

Then, in Sec. IV D 2 we focus on the so-called Toulouse 
point a — 1/2 of the ohmic model. For a single spin, one 
can solve for the dynamics exactly and ((T^(i)) exhibits 
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FIG. 13: Comparison of spin dynamics between the ohmic sin- 
gle spin-boson (thin lines) and ohmic two-spin boson model 
with Kr = (symbols). Different curves correspond to dif- 
ferent values of dissipation strength a. The Ising interaction 
is chosen accordingly to he K — Aauic- Other parameters are 
Ai,2 = A = 0.1 ojc, £1,2 = e = 0. The upper panel shows the 
spin dynamics in the coherent regime < a < 1/2. Here, 
the two-spin oscillations have a slightly larger frequency and 
are stronger damped than the single spin oscillations (see also 
Fig. 15 ). The lower panel displays the dynamics for stronger 
dissipation a > 1/2. For 1/2 < a < 1 both systems display 
incoherent decay. For even larger values of a, we observe that, 
in contrast to the single spin-boson model which localizes at 
a — 1, the two-spin dynamics remains incoherent at least up 
to a = 1.5. This is in agreement with the phase diagram of 
Fig. 2. 



pure exponential decay. The exact solution can most 
easily be derived by employing a bosonization mapping 
to a non-interacting fermionic resonant level model. In 
the two-spin case, we explicitly show that this mapping 
does not lead to a non-interacting fermionic model, which 
hence cannot be solved exactly. Further, our numerical 
results prove that the two-spin dynamics at a = 1/2 
differs slightly from the single spin case. We associate 
this with the influence of the retarded part of the bath 
induced Ising interaction which is still present even at 
Kr^O [see also Eq. (10)]. 



1. Weak dissipation: Quality factors and the 
Non-Interacting Blip Approximation (NIBA) 

In this Section, we derive a decoupling approximation 
that allows us to qualitatively understand the spin dy- 
namics for Kr = 0, but not necessarily small a. In 
this approximation, the two spins decouple completely 
for Kr — 0, and their Heisenberg equations of motion 
are identical to the ones of the single spin-boson model 
in the well-known NIBA.''''' 

Our starting point to investigate the dynamics is the 
polaron transformed Hamiltonian in Eq. (4). For zero 
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FIG. 14; Comparison of spin dynamics between the subohmic 
single spin (thin lines) and two-spin boson model with Kr = 
(symbols). Different curves correspond to different values of 
dissipation strength a. Ising interaction is chosen accordingly 
to be A" = SauJc- Other parameters are Ai,2 = A = 0.1 luc, 
£1,2 — e = 0. The upper panel shows the spin dynamics for 
weak dissipation up to a = 0.05. The two-spin oscillations 
have slightly larger frequency and are stronger damped com- 
pared to the single spin results. The lower panel shows the 
case of strong dissipation a > 0.05. The single spin-boson 
model localizes at Uc = 0.107. In contrast, the two-spin bo- 
son model remains delocalized along Kr = 0, which is in 
agreement with the phase diagram of Fig. 3. 



detuning and at Kr = 0, it reduces to 

2 

2 



^ = E^«^'''+h-c-) + E^4^'=' 



(40) 



i=i 



fc>o 



where fi — —i^^. —{bl — bk)- The Heisenberg equation 
of motion for cr^(t) with j = 1,2 reads 

a] (t) = -i A,a+ (t)e*^^(*) -t- h.c. . (41) 

It contains the elements a,, (t) which are given by 



4W 



iA,- 



dsa|(s)e-*"(^) 



(42) 



and a- ~ (o't)*- Inserting Eq. (42) in Eq. (41) yields 



A2 /■* 

^Ut) = -^ / ds[a^(s)e*"(*)e-*^^(") + h.. 
2 J_^ 



(43) 



Note that the two spins are still coupled to each other via 
the time-dependent bath operator fl{t). This coupling 
describes the retarded part of the bath induced Ising in- 
teraction that we have seen already in the spin effective 
action of Eq. (10). If we neglect this interaction, the two 
spins decouple from each other. 

More formally, we employ two approximations which 
are known to be equivalent to the NIBA in the single 
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FIG. 15: Quality factor of oscillations in the ohmic (upper 
panel) and subohmic (lower panel) single and two-spin bo- 
son models. TD-NRG results for the single (two-spin) boson 
model are shown as crosses (open squares). Quality factors 
derived from the Bloch-Redfield approach (see Sec. IV B) are 
displayed as filled squares. The solid line denotes the quality 
factor derived from the NIBA approximation of Eq. (45). 



spin case.'° First, we assume that the bath evolves freely 
bk{t) = fofc(0)e~*"'=*, neglecting any backaction of the 
bath on the spins. The reduced density matrix of the 
bath remains unaffected by the spins. Second, we trace 
out the bath degrees of freedom in a weak-coupling sense 

TrB[e^"«e-^"(^)] = exp{^ [*Qi(i - s) - Q2it - s)] } , 

(44) 
which contains the bath correlation functions Qi{t) = 
L duj J (uj) sin Ldt and Q2{t) — L dujJ{Lo)[l — cosujt]. 
As a result, the two spins are now completely decoupled 
from each other and their dynamics is described by 



^l(i) 



-A? 



ds 



■(s) cos 



Qi{t-s) 



-'Q2{t-S)/TT 



(45) 
Eq. (45) is known to describe the dynamics of the single 
spin-boson model in the famous NIBA.^^'" It can readily 
be solved by Laplace transformation. We refer to Refs. 2 
and 3 for details. 

In the ohmic case, we find from Eq. (45) that 
the spin undergoes damped oscillations for < 
a < 1/2. The frequency of the oscillations is 
given by wniba = Acff cos 2(fr^ ^'^d the damping 
rate reads Fniba = Acff sin 2(f°a) ' where A^s = 
[r(l - 2a) cos7ra]i/2(i-")A(A/tJc)"/(i-") is a renormal- 
ized tunneling element. The quality factor of the damped 



oscillations thus reads WNiBA/rNiBA 



cot 



2(l-c 



The 



NIBA also predicts an incoherent contribution to the 
spin dynamics, which is absent in the TD-NRG results. 
At a = 1/2, Eq. (45) predicts purely exponential re- 
laxation ((T^(/;)) = exp[— Fi] with a decay rate given by 

F = Aeff(a = 1/2) = f^. We refer to Ref. 9 for a de- 
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tailed analysis of the single spin-boson dynamics within 
the TD-NRG. 

In Fig. 15, we present a comparison of the quality 
factor of the oscillations w/F for the single and two- 
spin boson systems as computed by TD-NRG and Bloch- 
Redfield. It is obtained by fitting the numerical results to 



the function erf 2(i) 



-rt 



cos{ujt). In the ohmic case, we 



also include the prediction from the NIBA, which agrees 
with the TD-NRG results of the single spin-boson model. 
In general, we observe that the quality factor is smaller 
for the two-spin system. The Bloch-Redfield approach 
yields accurate results only in the ohmic case for small 
a. It fails completely in the subohmic case due to the 
increased spectral weight of slow oscillator modes, even 
for weak dissipation. 

If we increase the dissipation strength further to values 
a > 1/2, we observe an important difference between the 
single and the two-spin boson models. The two-spin bo- 
son model does not enter a localized phase for increasing 
values of a along the line K,. = 0. 

For an ohmic bath, spin transitions occur even for a 
dissipation constant as large as a > 1.5 (see Fig. 13). 
This is in stark contrast to the single ohmic spin-boson 
model, which displays a localization phase transition at 
a critical dissipation strength of ac = I + 0{ — ). This 
explicitly shows that the approximations that lead to 
Eq. (45) even fail to give the correct qualitative dynam- 
ical behavior for stronger coupling a. 

In the subohmic case, the NIBA cannot be justified 
and erroneously yields localization for all a > O," while 
TD-NRG results for the single spin-boson model show 
that the system remains delocalized up to a finite critical 
value of a.''" Again, we find in Fig. 14 that unlike the 
single spin-boson system, which localizes at a value of 
ac = 0.107 (for A/uJc = 0.1 and s = 1/2), the two-spin 
boson model always remains delocalized for Kr = 0. 

We finally want to emphasize that the NIBA in the 
single spin-boson model breaks down for finite bias e.'^ 
Thus, it cannot be applied to the two-spin problem away 
from the line Kr = 0, since the Ising interaction acts as a 
mutual bias between the spins. One common approach is 
to account for the interblip correlations up to first order 
in the spin-bath coupling a.'^''"*^ This procedure, however, 
is equivalent to the perturbative Redfield approach that 
we have discussed in Sec. IV B. 



2. Toulouse point: relation to the single spin case 

In this Section, we investigate the dynamics of the 
ohmic two-spin boson model at the special parameter 
point Kr = and a = 1/2. 

It is well-known that one can map the ohmic single 
spin-boson Hamiltonian in the scaling limit A/ujc ^ 1 
onto a fermionic resonant level model using bosonization 
and refermionization techniques.^ The fermionic model 
describes a localized level (dot) that is coupled via tun- 
neling to a lead of free spinless fermions. In general the 
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FIG. 16: Exponential decay of {a^{t)) at the Toulouse point 
a = 1/2 in the single spin-boson model. Symbols are TD- 
NRG results for different tunneling amplitudes A/ojc which 
result in different decay rates F = ^^- . Solid lines are fit with 
fi{t) = exp[— a^rt] using fit parameters a^ given in Table II. 



resulting model contains a Coulomb interaction term be- 
tween the fermions on the dot and the ones in the lead. 
At the special (Toulouse) point of a = 1/2, however, 
this interaction vanishes, and the fermionic model can be 
solved exactly, also in noncquilibrium situations. ^'^'^■^'^'' 
For a spin that is initially polarized along the z- 
direction, one finds purely exponential relaxation for 
t > (and e = 0)": 



{a^{t))=cxp[-Tt] 



(46) 



with decay rate F = 7rA^/2wc. It is worth noting that the 
NIBA predicts the same behavior, since it becomes exact 
at the Toulouse point of the single spin-boson model. ^ 

To prove the validity of the TD-NRG method in this 
strong coupling regime, we compare in Fig. 16 our nu- 
merical results for a single spin-boson model with the 
exact solution of Eq. (46). We observe that the decay is 
indeed purely exponential, and the decay rate is given by 
F = 1^ in the scaling limit A/ujc -^ (see Table II). 

From Eq. (45), we expect a similar behavior for the 
two-spin boson model at the (generalized) Toulouse point 
Kr — and a ~ 1/2. Indeed, as shown in Fig. 17, we 
observe that (17^2(0) decay purely exponentially in the 
two-spin case as well. The decay rates of single and two- 
spin models, however, are slightly different. We find in 
Table II that the decay rate of the two-spin boson model 
is about twice as large as the decay rate for the single 
spin-boson system. The difference of the decay rates is, 
again, due to the retarded part of the bath induced Ising 
interaction neglected in the derivation to Eq. (45). We 
will qualitatively explain the factor two difference below, 
using a mapping to a fermionic resonant level model. 

One might ask whether the two-spin boson model can 
also be solved exactly via the bosonization mapping to a 
fermionic resonant level model. For two spins, however. 
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FIG. 17: Exponential decay of (o-f 2(i)) at the generalized 
Toulouse point K = 2a;c, Q = 1/2 in the two-spin boson 
model. Symbols are TD-NRG results for different tunnel- 
ing amplitudes Ai,2/^c, which result in different decay rates 
ri,2 = ttAi 2/(2cJc)- Since Ai — A2, we observe Fi = F2. 
Solid lines are fit with fi{t) — exp[— 6iFt] using fit parameters 
bi given in Table II. 



it turns out that the fermionic model remains interacting 
at the Toulouse point, and thus cannot be solved exactly. 
As we show in detail in the Appendix, the (additional) in- 
teraction term is proportional to the tunneling elements 
Ai^2j which describe tunneling between dot and lead in 
the fermionic model. Since we are interested in a solu- 
tion that is non-perturbative in Ai 2 we cannot treat this 
additional term as a weak perturbation. 

Specifically, the ohmic two-spin boson Hamiltonian in 
Eq. (1) can be mapped to a fermionic resonant level 
model with two energy levels on the dot. The mapping 
becomes exact in the scaling limit Ai^2/<^c — )■ 0.^ As we 
derive in the Appendix, the resulting fermionic model is 
described by the Hamiltonian 



fc>0 j = l 

-V[{1~ «)ni4V'(0) + (1 + i)n2d\i^{0) + h.c] 
2 2 



j=i 



i=i 



KRL{d\d,-^)[dld2-l). 



(47) 



Here, Ck annihilates a spinless fermion of momentum k 
and energy cok = vpk in the lead {vp is the Fermi veloc- 
ity), and one defines ip{0) — L^^^^ ^j, c^, where L is the 
length of the lead. The colons denote normal ordering 
:^t(o)^(0):= i/jt (0)^(0) - (?/'^(0)V'(0)). The spin opera- 
tors have been expressed in terms of fermionic operators 
on the dot using the renowned Jordan- Wigner transfor- 
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TABLE II: Fit parameters for single spin-boson model {ai} 
and two-spin boson model {bi} with standard error a{ai) and 
a(6i), for different values of A/luc, or F = tvA'^ /2ujc. We fit 
the TD-NRG results to an exponential decay function which 
reads fi(t) = exp[— a^Ft] for the single spin-boson model, and 
fi{t) = exp[— feiFt] for the two-spin boson model. 



mation in a symmetric form' 



6.3,75,76 



"2 



= [1 - (1 - i)n2]di 

= [I - {1 + i)ni]d2 



(48) 
(49) 
(50) 



where rij = dljdj are the dot occupation number opera- 
tors. The parameters in Hj^l can be expressed in terms 
of the spin-boson parameters as 



2TrpV 



^A2 



pU 



1 



2a 



2wc 2 

Krl=K + 2uJc{1-2V2^), 



(51) 



where the fermionic density of states is defined as p = 
l/{2TrvF)- The bias field ej of the spin-boson model cor- 
responds to the energy of the dot level j with respect to 
the Fermi energy of the lead. 

The last two interaction terms vanish at the Toulouse 
point: U = Krl = for a = 1/2 and AT = 2ujc- The 
term, V [(1 - i)ni4V'(0) + (1 + i)n2d\ip{0) + h.c] , how- 
ever, is proportional to the dot-lead tunneling and thus 
remains. It arises due to the Jordan- Wigner string that 
accounts for the distinct commutation rules of fermions 
and spins at different sites. 

The dynamics of Eq. (47) cannot be solved exactly. 
Nevertheless, we can use the fermionic description to 
qualitatively understand that the decay rate of the two- 
spin boson model is about twice as large as in the 
single spin case. To this end, we introduce the sym- 
metric and antisymmetric combination of dot operators 
D, = [di + d2]/V2 and £»„ = [di - daj/v^. The 
occupation numbers can then be expressed as ni^2 = 
i [DlD, + DlDa± DlDa ± DlD,], where the upper' sign 
refers to ni. At the Toulouse point, the Hamiltonian then 
takes the form 



i?RL = Ho + E{ns + na) + AE{DlDa + h.c.) 
+ V2V{{Dl - Dlua - ^-Dln,)^(0) + h.c.} 



(52) 
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which contains the energy sum E = (ei + e2)/2 and 
difference AE — (ei — 62) /2. We write n^ = DlDg, 
Ua = D\Da and denote the free part of the lead electrons 
as Hq = 'vp'Y^^yQkc\ck- Both symmetric and antisym- 
metric state have the same energy E^ and the original 
energy level difference translates into an effective tunnel- 
ing coupling between them. 

For AE = and initially empty dots such that Ug = 
Ua = 0, the antisymmetric state decouples from the sys- 
tem completely. The tunneling coupling between the 
symmetric state Dg and the lead, however, is stronger 
than for each individual level di2- It is given by \/2V 
instead of V [see Eq. (47)], and the level will therefore 
fill twice as fast because T ^ V^. As soon as Ug > 0, the 
antisymmetric state couples to the lead as well, and in 
equilibrium one finds that (rij) = (ria) — 1/2 for E = 0. 
For AE — 0, symmetry requires that (ni) — (712) and the 
expectation values {D\Da) and {D\Ds) are thus purely 
imaginary. It then follows that (ni) = (712) = 1/2 and 
(crj) ~ (cr|) = in equilibrium. For AE 7^ the level 
correlations acquire a finite real part which gives rise to a 
difference in the level occupations (711,2) in equilibrium. 
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FIG. 18: Spin dynamics (erf 2(i)) for different values of a 
in the regime of strong spin-bath coupling for ohmic (upper 
part) and subohmic bath with s = 1/2 (lower part). Other 
parameters read A = O.lajc, K — 0.2uJc, £1,2 = 0. For this 
choice of K the localization phase transition occurs at Uc ~ 
0.25(0.043) in the ohmic (subohmic) system. 



E. Strong spin bath coupling 

In this Section, we focus on the regime of strong spin- 
bath coupling, where perturbative approaches are not 
applicable. We thus use the TD-NRG to calculate the 
spin dynamics, and focus on the differences between the 
case of an ohmic and a subohmic bath. We find that, 
qualitatively, the behavior in the two-spin boson systems 
resemble the one known from the (respective) single spin- 
boson model. For an ohmic bath, we observe in the upper 
part of Fig. 18 that the coherence of oscillations is lost 
above a certain bath coupling strength, roughly given by 
ac/2, where adK, Ai^2) denotes the critical value above 
which spin transitions are completely suppressed (local- 
ized regime of the phase diagram in Fig. 2). 

The situation is completely different for a subohmic 
bath as shown in the lower part of Fig. 18. Here, oscil- 
lations persist even into the localized region. This phe- 
nomenon was only recently discovered in the single spin- 
boson model '''"'^' and we confirm that it also holds in the 
two-spin case. 

This qualitative difference between the ohmic and sub- 
ohmic models at which point the coherence of the spin 
oscillations is lost (as a function of a and K), is also re- 
flected in the behavior of the static entanglement entropy 
(see Sec. Ill A). 



F. Generation of highly entangled steady states 

In this Section, wc show that the two-spin boson model 
displays interesting steady states for certain initial prepa- 
rations. In this state, the spins are entangled with the 



bath while maintaining coherence between different spin 
configurations. 

Let us ask the question what happens if we polarize the 
spins initially in an antiferromagnetic configuration such 
as Iti) in a region of the phase diagram where the ground 
state phase is localized. At i = 0, we then turn off the 
external bias fields completely, and follow the evolution of 
the spin reduced density matrix ps over time. Note that 
the system can only localize in one of the ferromagnetic 
spin states {| tt))] ii)} as discussed in Sec. II B. We 
calculate ps{t) using TD-NRG and observe that after a 
time of the order l/F = 2uJc/{t^A'^) the system reaches a 
steady-state where the spin reduced density matrix takes 
the form 



/l o\ 



Ps,^ 



1 





1 -1 


4 





-1 1 




\o 


V 



(53) 



where we use the standard basis {|tt), Iti), lit), III)}- 
With a probability of j, the spins are thus localized in 
one of the ferromagnetic spin states {|tt), lii)}; and with 
a probability of | the spins are in the spin singlet state. 
The entanglement entropy £, which is a measure of the 
entanglement between spins and bath, is nonzero in this 
state. Specifically, S{ps.ss) — | from Eq. (53). 

We can easily understand this form of the steady state 
by writing the initial state in terms of the singlet state 
|5 = 0,m =: 0) == [|tl> - |it)]/\/2 and the triplet state 
\S^l,m^Q) = [\n) + \m/V2as 

in) = ^{\S = l,m = 0) + \S^O,m^ 0)) . (54) 
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Whereas the singlet state does not couple to the bath 
at all, the triplet state localizes in one of the two ferro- 
magnetic configurations. In this steady state, the spins 
are highly entangled with the bath modes, while devel- 
oping and maintaining coherence between the two anti- 
ferromagnetic spin configurations. 



V. CONCLUSIONS 

We have presented an extensive study of a system of 
two Ising-coupled quantum spins in contact with a com- 
mon bosonic bath. We have investigated several distinct 
equilibrium and nonequilibrium situations, both for the 
case of an ohmic as well as a subohmic bath. Employ- 
ing the bosonic numerical renormalization group (NRG) 
and its recently developed time-dependent version (TD- 
NRG), we were able to describe the complete range of 
parameter space, from weak-to-strong coupling. We have 
applied a variety of different analytical approaches to 
comprehend, interpret and validate the numerical results. 

Using NRG we have calculated the ground state phase 
diagram of the model for s = 1/2 and s = 1. We find 
a striking asymmetry in the behavior for ferromagnetic 
{K < 0) and antiferromagnetic {K > 0) Ising coupling, 
which we have understood as being the result of the fact 
that the system only localizes in a ferromagnetic spin 
configuration {|tt), |ii)}- 

Let us briefly comment on the case of an SU(2)- 
symmetric spin-spin interaction ^cri ■ c^- First, due 
to the fact that the spin couples to the bath via its 
cr^-component, only the Ising component of the spin- 
spin interaction becomes renormalized by the bath K^. 
This generates an anisotropic XXZ-coupling of the form 
^(crfcrf -(- a\al) + ^crfcr|, where K^ = K. To ar- 
gue that the physical properties of a such a model are 
quite distinct from the Ising case where K^ = 0, we 
employ the well-established mapping to a two-impurity 
Kondo model. ■^''''^' In our case, it turns out that the trans- 
verse part of the coupling is invariant under this map- 
ping, and the Ising component renormalizes to K^ — 
K + 4lJc(1 - 2yja). As shown in Refs. 36,78,79 for the 
isotropic two-impurity Kondo model, the behavior of a 
Kondo system with K^ ^ greatly differs from the 
pure Ising case. In particular, in the absence of particle- 
hole symmetry the phase transition is replaced by a 
smooth crossover, whereas in the presence of particle- 
hole symmetry a phase transition occurs, but it is not of 
Kosterlitz-Thouless type.^'^ A detailed (numerical) anal- 
ysis of the SU(2)-symmetric two-spin boson system is left 
to further studies. 

Here, we have then investigated the behavior of the 
Ising two-spin boson system close to the localization 
phase transition, which is in different universality classes 
for s = 1 and s < 1. In the ohmic case, we find that 
coherence in the ground state is lost prior to localiza- 
tion. This is reflected in a plateau in the entanglement 
entropy, which describes the entanglement between spins 



and bath. Eventually at a critical coupling strength, the 
spin is localized where the entanglement entropy quickly 
drops to zero. We have reported that the size of the 
plateau shrinks considerably for larger values of the Ising 
coupling constant K > Wc, indicating that, in this case, 
spin coherence is lost only close to the phase transition. 
Whereas the transition is in the Kosterlitz-Thouless uni- 
versality class for the ohmic system, it is of continuous 
type for a subohmic bath, where we have studied the 
scaling of the spin magnetization {a^ 2) close to the tran- 
sition. We have extracted critical exponents using NRG 
and compared them to analytical mean-field exponents. 
The agreement is reasonable though not perfect, which 
shows that NRG goes beyond the mean-field approxima- 
tion that we have used. 

In the last part, we have discussed a number of differ- 
ent nonequilibrium scenarios. First, we have investigated 
the exactly solvable case of zero transverse fields where 
TD-NRG agrees perfectly with the exact analytical so- 
lution. For weak-spin bath coupling, we have provided 
quantitative limits on the applicability of the commonly 
used perturbative Bloch-Redfield method. 

The coupling to the bath can be exploited to dynami- 
cally synchronize spin oscillations, which can prove useful 
in cases where a direct coupling between the spins is un- 
available. Since the bath induced Ising coupling scales 
with the large bath cutoff frequency Wc, synchroniza- 
tion even occurs at small a where decoherence is weak. 
Nevertheless, this phenomenon could not be observed 
within the perturbative and Markovian Bloch-Redfield 
approach. 

We have then investigated the dynamics of the two- 
spin boson model for K^ — 0, and have pointed out simi- 
larities and differences with the case of a single spin. We 
have derived the mapping of the two-spin boson model 
to a fermionic resonant level model, which contains two 
levels on the dot. In contrast to the single spin case, 
this model remains interacting at the Toulouse point due 
to an additional interaction term that arises from the 
Jordan- Wigner transformation of the spins. 

We have further studied {af 2(0) i'^ the crossover from 
weak to strong coupling where perturbative approaches 
cannot be applied. For strong coupling we have found 
that while spin transitions do not occur in the localized 
regime for the ohmic system, coherent spin oscillations 
persist into the localized regime for a subohmic bath. 

Finally, we have shown that the system features an 
interesting steady-state if we initially prepare it in an 
antiferromagnetic spin configuration within the localized 
regime. In this state, the spins are highly entangled with 
the bath degrees of freedom, and still develop and main- 
tain coherence between different spin states. 
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Evaluating this transformation term by term, and per- 
forming the sum over wavevectors 



L 2^ 



1^ = ^, one finally obtains the Hamiltonian 



fc>o' 



-ak 



Appendix: Mapping of the two-spin boson model to 
the fermionic resonant level model 



In the Appendix we provide details of the mapping of 
the two-spin boson model to the fermionic resonant level 
model. Due to the Jordan- Wigner string, the correspond- 
ing fermionic model remains interacting at the Toulouse 
point in the case of two spins. 



2,36,75 



Using bosonization techniques,^'''"'"' one can map 
the two-spin boson Hamiltonian in Eq. (1) with an 
ohmic spectral density J{uj) — 27raa;exp[— w/wj onto a 
fermionic resonant level model, which describes a central 
region (dot) coupled via tunneling to free spinless elec- 
trons in the lead. The number of spins in the spin-boson 
model is equal to the number of levels on the dot, and 
the number of bosonic baths is equal to the number of 
leads. 

Our starting point is the two-spin boson Hamiltonian 
in Eq. (1) 



H: 



SB 



y ' 2 J' ^ 2 ^ 

3 = 1 


+ 


fc>0 


M 


+ bk) 


+ ^^l'^2 


+ y ^ ujkblbk . 

k>0 










(A.l) 



To obtain the mapping to the resonant level model (and 
similarly to the Kondo model), where the bath consists 
of free fermions, we choose the oscillator dispersion to 
be linear ujk = vpk, with Fermi velocity vp, and the 
coupling constants 



Xk 



'a2vF 



Trk 



1/2 



'fc/2ci 



(A.2) 



The bath spectral density J(oj) = '?rX]fc>o '^t^i'^ — ^k) is 
then of ohmic form J{lo) = 27raa;exp(— w/ojc) up to an 
exponential cutoff at ujc- If we insert this into Eq. (A.l), 
the spin-bath coupling term becomes 



E-lEl 



2avp 



i=i 



fe>0 



2ttL 



1/2 



-"""/^{bl + bk), (A.3) 



where we have defined the small distance cutoff a = 
We now apply a unitary (Luther-Emery) transforma- 



tion to the Hamiltonian: -ffsB = U^H^^U ^ where 



t/-y = exp 



E=i<? with 



x^rA, 



i?SB = ^F E k^l^^ + E{l"«^'' + ^■^■) + i-3 



k>0 



J = l 



nvp 



(727 - V2a)a] ^ e'^ 



-ak/2 



k>0 



2ttL 



1/2 



(bk + bl)} 



+ {K + AlOc7^ - SucVaj) -^ 



(A.5) 



One can show that a particular combination of the Bose 
operators bk , b^ can be made into an anticommuting 



Fermi field tp{x) 



: expj{x) with 



jW = E' 

k>0 



-ak/2 



r27r 



kL 



1/2, 



{bke 



ikx 



t ^ — ikx\ 



Ke 



(A.6) 



The coefficients have been chosen such that 
[j{x),j{y)] — —iTTsign{x — y) for a — > and thus 
{exp±j(x),exp±j(2/)} = 0, for x ^ y. Choosing 
7 = l/v2, one can thus identify the exponential 
exp[^/-\/2] which multiplies <7^2 ™ Eq. (A.5) as a local- 
ized spinless fermionic field ^(0) — (27ra)~^'^ exp^/\/2. 
The bosonic oscillator degrees of freedom are then 
interpreted as the density excitations p{k) = X^pcLfc^pj 



p{—k) ~ p^{k) of the fermions ^Ij{x) = L ^^^^ 



271- 
kL 



1/2 



fc>0 CfcC 



p{-k). 



via the bosonization identity bk ~ 

Using refermionization we can replace the free bosonic 
with a free fermionic Hamiltonian vp X]fc>o ^^I^fe ~^ 
VF Y.k>o ^4cfc and 



^ e-"'=/2 [ 



fc>0 



k 
2^ 



1/2, 



{bk + b^ 



J2^-^[p{-k)+p{k)]^:i;\0)i;{0):, (A.7) 



k>0 



where : V^(0)V'(0) := ipHO)ip{0) - {^pHO)ip{0)) denotes 
normal ordering. Finally, we write the spin operators 
in terms of fermionic dot operators using the Jordan- 
Wigner transformation (in symmetric form) 



T-^ = [1 - (1 -i)n2]di 


(A.8) 


^2^ = [1 - {l + i)ni]d2 


(A.9) 


a|-2nj-l,forj = l,2. 


(A.IO) 



E^- 

fc>0 



ak/2 



in 
kL 



1/2, 



(bk 



(A. 4) We note that a less symmetric form of the transformation 
is equivalent. The Hamiltonian (A.l) thus reads in terms 
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of fcrmionic fields as 

2 

fe>o i=i 

- [^2(1 - i)ni4V'(0) + Vi(l + i)n2d\tl^{0) + h.c. 



KRL{d\d^-^)(dld2-l), (A.ll) 



with tunneling coupling constant Vj 



2 



1/2 



dot- 



'J ^ ~2~ l^ptlJI 

lead interaction [/ = (1 — \/2a)/2p and onsite coupling 
-K^RL = K+2ujc{l—2y/2a). The interaction parameters C/ 
and K^i^ vanish at the Toulouse point a = 1/2 and K = 
2ujc- The additional interaction term [V2{l—i)nid*2tp{0) + 

Fi(l + i)n2d\tp{0) + h.c], however, is proportional to 
the tunneling couplings Vj and remains present at the 
Toulouse point. As a result, the fermionic model cannot 
be solved exactly. 
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